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Abstract 

In this work we combine information from relic abundance, direct detection, cosmic microwave 
background, positron fraction, gamma rays, and colliders to explore the existing constraints on 
couplings between Dark Matter and Standard Model constituents when no underlying model or 
correlation is assumed. For definiteness, we include independent vector-like effective interactions 
for each Standard Model fermion. Our results show that low Dark Matter masses below 20 GeV 
are disfavoured at the Scr level with respect to higher masses, due to the tension between the relic 
abundance requirement and upper constraints on the Dark Matter couplings. Furthermore, large 
couplings are typically only allowed in combinations which avoid effective couplings to the nuclei 
used in direct detection experiments. 
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I. INTRODUCTION 


The search for Dark Matter (DM) in the form of thermal relics represents one of the most 
active lines of research in astro-particle and particle physics. Indeed, there is an overwhelm¬ 
ing number of dedicated experimental searches for DM, most of them concentrating on the 
so-called Weakly Interacting Massive Particles (WIMPs) paradigm [1]. These are classi- 
hed into three different categories: (1) indirect detection searches, where the DM would 
annihilate or decay into SM particles which can be detected, (2) direct detection searches, 
where the DM would scatter against the protons and neutrons in a detector, producing an 
observable recoil, and (3) collider searches, where the DM would be produced in high-energy 
collisions, thus leading to events with missing momentum.^ 

In light of all these searches it is essential to look at the global picture of where we 
stand concerning the WIMP scenario. This is certainly a rather complicated task, given 
the enormous zoo of models available in the literature with an associated plethora of free 
parameters. Even for a single model, brute-force scans of the corresponding parameter space 
represent a signihcant computational challenge. Consequently, the usual practice is to either 
rely on simplified realisations of those models or to £x some of the free parameters, thus 
scanning only hyper-surfaces in the model parameter space. ^ 

On the other hand, Monte Carlo methods constitute an efficient alternative to scan over 
a multi-dimensional parameter space. However, for WIMPs, most of the Monte Carlo scans 
are performed in the context of supersymmetric models, where the neutralino is usually 
selected as the preferred DM candidate (e.g.. Refs. [IMH, however see na as an example 
of a non-SUSY search). 

The main goal of this work is to explore the present status of our knowledge of DM cou¬ 
plings to the different SM constituents in as much generality as possible. In order to gain in 
model independence, it is thus interesting to be more agnostic about the DM interactions 
with the SM constituents and parametrize them through an effective field theory (EFT) 
approach. In particular, this approach can reveal how constrained the DM interactions with 

^ A complementary search strategy is the search for DM self-interactions which would impact structure 
formation as well as stellar evolution in particular scenarios (see e.g. EHS]). Since the focus of this work is 
to probe the interactions between DM and visible particles these constraints will not be considered here. 
^ The literature in this respect is quite vast. Some examples relevant to our work are EHU, in the context 
of effective field theory (see below). 
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the different SM fields are in the light of currently available data, when no further assump¬ 
tions are made. Since allowing all type of Lorentz structures for these effective operators 
would provide too much freedom with hundreds of operators that cannot be bounded, we 
rather choose to exemplify the different interactions between a Dirac fermion DM and the 
SM by operators of the form 

Ci,p ixi^x) , ( 1 ) 

with independent coefficients for all SM particle species fi with chiralities P = Pp,Pp. 
This type of interaction is well motivated by scenarios where the DM communicates with 
the SM via a vector portal. Moreover, scalar or tensor interactions typically require Higgs 
insertions, leading to higher dimensional and hence more suppressed operators. Addi¬ 
tional dimension six operators with a vector-like structure involving the Higgs, such as 

^ _ y 

iXX^xH'^ D could be included as well. For DM heavier than the Higgs, this opera¬ 
tor would provide an additional annihilation channel relevant for indirect searches and relic 
abundance constraints, via xx hh. However, as we will see, above 100 GeV the global up¬ 
per bounds for any individual coupling correspond to the ones obtained from relic abundance 
alone (with no constraint for the top quark coupling, since that channel is not kinematically 
allowed). Therefore, we expect a similar behaviour for the Higgs, that is, no constraints for 
DM masses below the Higgs mass and recovering the relic abundance constraint for higher 
DM masses. Even with this restriction, if independent couplings for all SM constituents are 
allowed, 15 independent parameters remain to be constrained. Thus, given the large dimen¬ 
sionality of the parameter space, we make use of a nested sampling Monte Carlo algorithm 
to scan it efficiently. 

When constraining the DM EFT parameter space, we consider bounds from all types of 
experiments where a WIMP signal is being actively sought for, i.e., direct detection (namely, 
LUX |T6] and EDELWEISS [Hj), indirect detection (AMS [18] positron fraction data and 
Fermi-LAT data for dwarf galaxies mi), cosmic microwave background (CMB) and relic 
density constraints from Planck [20], and monojet and monophoton searches in colliders 
(from LHC [21] and LEP data [H])-^ 

^ Contrary to the cases of direct and indirect detection strategies, for collider searches the use of the EFT 
framework -mainly at the LHC- may not be optimal. We take this issue into account when recasting 
the limits coming from LHC. See Sec |HI| Also, note that other searches different from monojets can be 
interesting in the context of specific models (see e.g. Ref. |22) . 
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Additional information on some annihilation channels such as rr, bb, or uu come from 
the non-observation of indirect neutrino signals from the Sun in SuperKamiokande [ 23 ] and 
IceCube [ 23 |- However, the bound is only relevant if apart from a large DM coupling to the 
hnal state particle, a sizeable coupling to hrst family quarks is also present. Since we are 
assuming all couplings to be independent, this bound will always be satished in the £t by 
choosing a set of parameters with one of the relevant couplings being very small. Thus, the 
inclusion of these experiments would not change our results. 

We derive bounds on the coefficients accompanying each effective operator as a function 
of the DM mass, as well as bounds on the DM mass itself. These constitute the most general 
constraints which can be derived assuming that DM interacts with the SM as in Eq. Q. 

The outline of the work is as follows. In Sec. [IT] we describe the set of effective operators 


that will be jointly analyzed, introducing the parameters to be constrained. Sec. |III| lists the 
set of experimental constraints considered, and how these have been implemented. Some 
details regarding the numerical tools employed in the £t are explained in Sec. |IV[ Finally, 
Sec. [V| summarizes our results and we present our conclusions in Sec. |V^ 


II. THE EFFECTIVE FIELD THEORY FRAMEWORK 

In order to explore how much freedom the present global data on DM allows for its 
couplings to the SM, we will dehne a series of working models characterized by a set of 
independent effective operators, described by the following effective Lagrangian 

= {xrx)jf- ( 2 ) 

where the effective current is given by 

^ 7/^ I ’ I + CQi 7 m I ’ I + X] CfRfR'y/.fR, (3) 

£=e,^,T y J i=l J f 

the index i denotes the quark generations such that Ui = u, U 2 = c, = t, di = d, 
d 2 = s, and d^ = b, while the sum over / runs over all right-handed (RH) SM fermions. The 
coefficients cx are the couplings of the operators in the effective Lagrangian. As the effective 
operators are of dimension six, these coefficients will have a mass dimension of minus two. 
In expr.(|^, y represents a Dirac fermion DM. 
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This set of operators provide enough freedom so as to parameterize DM interactions with 
different strength to the various SM particles while, at the same time, keeping the set of 
operators at a viable level. Indeed, allowing extra operators would not only imply a more 
challenging numerical analysis, but would also be rather uninformative since mostly any 
value could be htted and particular UV completions tend to have a much more limited set 
of free parameters. Furthermore, for simplicity and due to the generally stronger constraints, 
we will not consider flavour-changing operators between the SM fermion bilinears [25] . De¬ 
spite these restrictions, 15 different operators fall into this category, parametrizing the DM 
interactions with the 3 quark and lepton doublets as well as the 3 singlets of up-type and 
down-type quarks and charged leptons. Starting from this general setup, we will also dehne 
more restrictive working models that can exemplify other interesting scenarios such as lep- 
tophilic, leptophobic or flavour-blind setups, characterized by different correlations among 
the couplings. The models considered in this work are the following: 

1. General model: The first model makes no additional assumptions regarding the coef- 
hcients cx, which are all allowed and free to vary independently. This represents the 
least restrictive choice possible in our given EFT framework, and includes a total of 
15 coefficients in addition to the mass of the DM. 

2. Flavour-blind: In this model, all operators involving fermions with the same gauge 
quantum numbers, e.g., all left-handed (LH) quarks, are assumed to also have the same 
couplings to DM and therefore have the same coefficients in the effective lagrangian. 
We are then left with five different coupling parameters which are related to those of 
the general model as 

^Le ^Q1 ^Q2 ^Q3 — Cq, 

^uR GR GRi GR Cgfi GRi ^eR GR- ( 4 ) 

For the latter three groups of coefficients, we will generally use the notation for the 
coefficient of the hrst generation to refer to it within this model. 

3. Family-oriented: In this model, we make the assumption that the effective couplings 
with DM are equal for all particles belonging to the same generation. This can be 
considered a quite crude proxy to flavour theories were the successive SM families are 
characterized by hierarchical couplings or charges in order to explain the observed mass 
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hierarchy between them, following the philosophy of the Froggatt-Nielsen mechanism^ 
[2MI]. This leaves only three independent operator coefficients, 


C-Le ^eR ^Q1 ^uR ^dR — *^Q2 ^cR ^sR — ^2y 


CLt = CrR = Cq 3 = CtR = CbR = C3. 


(5) 


This is the model with the fewest number of free parameters which we will consider and, 
as such, the correlations among the couplings will allow to obtain strong constraints, 
particularly for the coupling Ci. 

4. Leptophobic: This model assumes that DM does not have significant interactions with 
any of the leptons. As a result, all of the coefficients associated to operators involving 
leptons are set to zero, i.e.. 


( 6 ) 


cie — Cm — 0 , 


with ^ = e, fi, T. At the same time, no restrictions are imposed on the operators 
involving quarks and the model therefore contains nine free parameters. 

5. Leptophilic: In direct contrast to the leptophobic model, this model instead considers 
the situation where the only relevant DM interactions with the visible sector are those 
involving leptons. In this model, we therefore set all of the coefficients for operators 
involving quark helds to zero, i.e.. 


(7) 


CQi — CuiR — CdiR — 0 , 


where i indicates the generation. In this situation, we are left with six independent 
coupling parameters for the DM interactions with leptons. 

We wish to stress that these models are only meant to be phenomenological tools to assess 
our present global knowledge of how well constrained the DM couplings to the SM fermions 
are when not assumed to be universal or related through any particular UV completion. We 
will thus allow all couplings to vary freely when fitting the present data without questioning 
the apparent naturalness (or lack thereof) of the preferred regions. In particular, as we will 
see in Sec. |V| most models prefer 2 cqi + Cur + CdR ~ 0 so as to avoid the very stringent 

^ See [ini for a realisation in the context of dark matter. 
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limits coming from direct detection searches. Thus, even if these cancellations may seem 
unnatural, we will not avoid them by adding artihcial “naturalness priors” to guide the 
Monte Carlo in any way, in the spirit of allowing the Monte Carlo to choose the points in 
parameter space which provide the best £t to the data. Symmetry arguments could perhaps 
be invoked in a particular UV completion when trying to reproduce the best £t found in the 
effective description in a natural way. 

Moreover, these working models are not intended to be self-consistent low energy de¬ 
scriptions of any particular UV completion. In fact, due to the limited amount of operators 
considered, we are not allowed to take our effective theory description beyond tree level 
processes. Indeed, radiative corrections could induce other operators with interesting phe¬ 
nomenology.® For example, in the leptophilic model the lepton legs could be closed in a loop 
and, through the emission of a virtual photon, induce a coupling to hrst generation quarks. 
Therefore, the strong bounds from direct searches would also apply to these models and 
put a constraint on the lepton couplings (see Ref. [M])- In a similar fashion, other signals 
at the LHC such as dijet or dileptons could be generated from loop-induced 4-SM-fermion 
operators. Furthermore, when moving from the flavour basis to the mass basis, the DM 
couplings to the SM doublets will induce flavour-changing operators that, when the DM legs 
are closed in a loop, could contribute to the oscillations of neutral kaons or other FCNC 
processes, for which stringent experimental constraints would apply. In full consistency, 
these operators should have been included from the start, since they are compatible with 
the particle content and symmetries of the theory and they are required to renormalize the 
divergences stemming from the loop contributions. However, when doing so, new unbounded 
coefficients are introduced and, in order to avoid these stringent limits, the £t will prefer the 
points in the parameter space where these new coefficients exactly cancel the loop-induced 
contributions. Thus we will simply not consider these loop-induced constraints in our list 
of observables, since in our approach they would not imply additional constraints on our 
parameter space unless particular relations between the couplings are invoked. 

® For recent works in an EFT framework including one loop processes, see e.g. Refs. [321133j 
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III. CONSIDERED CONSTRAINTS ON DM 


In our numerical analysis, the DM contributions to the different DM observables are 
computed with micrOmegas [3S]- These are then compared to the current experimental 
bounds. A chi-square function is then computed for each observable independently. The 
total is obtained as the sum of all separate contributions 




( 8 ) 


where i is an index which runs over the observables and c is a vector containing the coupling 
coefficients of the model being tested. 

The contribution to the chi-square from each experiment is computed assuming that the 
limit on the physical quantity being bounded (he., relic abundance, the thermally averaged 
cross section, etc) by the experimental collaborations behaves like a gaussian. The chi-square 
can be expressed as a function of the number of events as: 



and ANgys refers to the systematic uncertainty on the number of events. Notice that, given 
the unfortunate lack of a positive DM signal in any of the observables described below with 
the exception of the relic abundance constraint, in this expression it has been assumed that 
the observed number of events corresponds with the expected background. In some cases, 
small upward or downward fluctuations over that background will be present and thus small 
deviations with respect to the scaling derived below can take place. However, notice that 
we will normalize this scaling with the data presented by each experimental collaboration 
itself. Thus, our function will always reproduce correctly the result obtained by the 
collaboration at the confidence level at which it is reported (typically 90 or 95%) and the 
deviations should be small for the purpose of our exploration as long as CL close to these 
limits are investigated. 

Experimental collaborations present their results as bounds on a given quantity at a 
certain confidence level (CL). This can be translated into a certain number of events 




allowed at that CL. For a bound at 68% CL, for instance, we get®: 


X68% — 1 — 


J^exp 


( 10 ) 


'sys 


By taking the ratio between Eqs. ([^ and (10), we get the following expression for the 


chi-square: 


= 


( 


N^h \ Nexp ^ ^ 


\N--p) + Nbg + ANsys ■ 

The expression above can be further simplihed in the two following limits: 


sys 


( 11 ) 


Statistically-dominated: in this case, the number of background events and the sys¬ 
tematic errors can be neglected, and we can therefore write the chi-square contribution 
as: 

]\fth 

( 12 ) 


^ - 


JSfexp 


Background/Systematics-dominated: in this case, the signal events can be neglected 
with respect to the background and/or the systematic error, and the chi-square sim- 
plihes to: 


- 


/ jY*^ Y 

I / 


(13) 


By making use of Eqs. (12) or (13), all the hner details of the detector response cancel 


in the ratio. We will therefore use the experimental bounds on the different observables 
reported by the collaborations instead of the number of events (which are not always publicly 
available). 

All the experiments considered in this work fall in the background/systematics-dominated 


case, and therefore we will apply Eq. (|13|). The only exception to this will be the case of dwarf 

These are dark-matter-dominated 


galaxies studied by the Fermi-LAT collaboration 
objects and constitute extremely clean probes to test for dark matter interactions with SM 
particles. We have explicitly checked the behaviour of the chi-square for a subset of the 
observed dwarf galaxies, using publicly available data from the collaboration from Ref. [19] , 
and we hnd that the behaviour is in between the two extreme cases identihed above. We 


therefore opt for the most conservative approach in this case, which corresponds to Eq. (12), 


® When the experimental result is quoted at a different CL, it can be simply rescaled to 68% CL within a 
gaussian approximation. 


9 














when using the limit provided by the collaboration (which has been obtained using the full 
data sample). 

In the following, we describe the full set of observables sensitive to interactions between 
the DM and the SM fermions which have been included in our simulations. 


Relic abundance 


The abundance of DM in the Universe is well-known from the PLANCK measure¬ 
ments m- At present, its central value is 


= 0.1187, 

X 


(14) 


with an error of an = 0.0012 at Icr. We will assume that our DM fermion y constitutes all 
of the observed relic density and that it is produced in the early Universe through thermal 
freeze-out. The predicted relic density D^^(c,m^) is computed for each set of values for the 
model parameters, and the corresponding contribution to the chi-square function is given 
by 


xUc,m^) 


V / 


(15) 


In general, within the thermal freeze-out scenario, the relic abundance of DM is mainly 
governed by the thermally averaged total DM annihilation cross section 


(av) = '^Wi{av)i oc 

i 


where 



(16) 


Here, i runs over all fermion helds contributing to the annihilation cross section, and Wi is 
a weight associated to the dimension of the SM gauge group representation of the corre¬ 
sponding fermion held. As we will show in the results section, this combination is subject 
to very strong constraints and dictates the possible ranges in which the coefficients may he. 
In particular, due to the hnite and non-zero relic abundance and the assumption of thermal 
freeze-out, this implies that at least one of the coefficients must be non-zero and that none of 
them can be too large. It should be kept in mind that, alternatively, under the assumption 
of freeze-in the correct relic abundance could be generated via extremely small couplings to 
all SM fermions [36] . 
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Direct detection 


The contribution to the chi-square from direct detection experiments is computed assum¬ 
ing that the limit on the spin-independent cross section given by the experimental collab¬ 


orations behaves like a gaussian and is background-dominated (see Eq. (13)). Therefore, 
the contribution to the chi-square coming from each direct detection experiment can be 
computed as 


Xlnic,m^)=(^) , (17) 

where we have expressed the limit at the 68% CL (after rescaling from the 90% quoted 
by the collaborations). Note that, for = cr^^, the appropriate value of the chi-square 
function and the experimental limit are recovered. 

For we implement direct detection constraints from LUX [16] and, so as to have 
an independent target material, EDELWEISS-II [I^. Our setup generally leads to a spin- 
independent contribution (except when cqi = —= —CdR) which will be constrained 
by the LUX and EDELWEISS-II results. For our set of operators, the spin-dependent 
contribution cancels out m- The spin-independent dark-matter-nucleus coupling will 
relate to the quark couplings as 


cn — ^[3^ Cqi -|- (% -|- Z)cur + {2A — Z)cdR], 


(18) 


where A is the mass number and Z the atomic number of the target nucleus. With this 
dependence on the couplings, there is a possible degeneracy in the DM-quark couplings for 
which the cross-section vanishes. For Z ~ %/2, this degeneracy occurs for 2cQi + CuR + CdR — 
0. Thus, direct detection experiments will only allow large couplings to the hrst generation 
quarks if this particular relation is fulfilled. As we will show in Sec. |V] this is clearly seen 
from the numerical results of our simulations. 


Cosmic Microwave Background (CMB) 

DM annihilations result in an injection of power into the Intergalactic Medium (IGM) 
per unit volume equal to EH Eg 

f 8^) = , (19) 
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where z is the redshift, Q^^o(Pc,o) is the DM abundance (critical density) today, (av) is the 
thermally averaged annihilation cross section, and the statistical factor ( = 1/2 corresponds 
to DM being a Dirac particle. On the other hand, CMB probes such as the WMAP and 
Planck satellites can set limits on the deposited power into the IGM around the CMB epoch 
(z ~ 1000), which is related to the injection power as [391 - HT] . 


dE 

dtdV 


dep 




dE 

dtdV 


inj 


( 20 ) 


where the efficiency function fj depends on the DM annihilation channel, yy —)■ PjPj- For 
this reason, experiments usually quote their limits through the quantity 


Pann-5^/, 2^' > (21) 

j ^ 

where {(Tv)j is the thermally averaged partial annihilation cross section into particles pj. 
This quantity has been constrained by Planck+lensing data at the 95% CL [20], here we 
rescale that bound to the 68% CL obtaining: < 1.7 x 10“^® cm®s“^CeV“^ so as to 

build our y^ function . In order to implement these bounds we use the tabulated values of 
fj from Ref. [2S], to compute compare it with the experimental constraints. The 

contribution added to the total y^ is given by 

Xcmb(c5'^x) ^ • (22) 

xPann / 


Positron fraction 

The measurements from AMS02 on the positron fraction [TH] are used to derive upper 
bounds on DM annihilation cross sections. The contribution to electron and positron fluxes 
$e±,DM from DM annihilation are computed using micrOmegas. The contributions coming 
from astrophysical sources, on the other hand, are parameterized by 

$e-,bg(F;) = + C<,E-^\ <h,+,bg(E) = (23) 

as in Ref. 021 , where the last term in the positron flux represents a point source with a 
hard spectrum cut Eg, while the other terms model the diffuse background. Both DM 
and background fluxes are affected by the solar modulation, which can be explicitly taken 
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into account by computing the flux at the top of the atmosphere (©) under the force held 
approximation: 

= + (24) 

where refer to the parameters accounting for the modulation. The differential positron 
fraction on top of the atmosphere is then given by 
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(|)® <j)© 

^e+,bg ' ^e+,DM 


)© I ,|,© I ,J)© I ,|,© 

e+,bg “ ^e+,DM ^ e-,bg “ ^e-,DM 


(25) 


The function is built by binning the predicted positron fraction and comparing to the 
experimentally measured values, 


bins 


xiMs(c,"ix) = 


7j -^exp,j 


(26) 


where the uncertainties in each bin ctj-j are taken directly from Ref. [18] by adding the 
statistical and systematic errors in quadrature. 

We have checked that allowing the electron flux to vary does not affect the fit in any 
substantial way and have therefore fixed the parameters for this flux to the values which 
provide the best-fit to Fermi electron flux data |13|. In particular, the value of the solar 
modulation parameter for the electron flux which gives a best-fit to the Fermi data is found 
to be equal to zero. Moreover, the solar modulation parameter for positrons (p'^ was found to 
have a negligible impact on the fit, while creating some numerical degeneracies. Therefore, 
we choose to fix this parameter to zero as well during the fit, as described in detail in App. [A] 
However, it should be noted that the solar modulation only affects the low energy part of the 
spectrum and therefore our analysis of AMS02 is expected to be accurate for > 10 GeV. 

In order to take into account the lack of knowledge on the astrophysical backgrounds, we 


profile over the positron flux parameters in Eq. (23) in our simulations,' without imposing 
any priors on them. This is done for each set of parameters {c} independently, as explained 
in detail in App. 


^ The minimisation procedure was performed with the GSL libraries |44j . 
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Gamma rays from dwarfs 


Dwarf spheroidal satellite galaxies (dSphs), being DM-dominated objects, constitute very 
clean laboratories to test DM interactions. A search for 7 -rays coming from Milky Way 
dSphs has been performed by the Fermi-LAT collaboration [19]. We make use of the present 
90% CL bounds (rescaled to the 68 % CL) on the thermally averaged cross section for DM 
annihilation from Ref. [19], which are derived under the assumption of 100% DM 

annihilation into some specihc channel XX PjPj- However, in our EFT approach, DM may 
annihilate to all fermionic channels with different branching fractions Rj, which will depend 
on the set of couplings for a given model {c}, and each particular channel j will contribute 
to 7 -ray emission with a different weight, given its different subsequent decay chains. This 
weight is inversely proportional to the hnal experimental bound that can be derived for that 
particular channel We may therefore recast the experimental bound on the total 

cross section as 




B, 


-1 


(27) 


As already mentioned, the contribution to the chi-square coming from dwarf galaxy mea¬ 
surements is not completely dominated by either background/systematics nor by statistics. 
Therefore, we conservatively assume that it is the latter which dominates the measurement 
in this case. The contribution to the is thus given by 


XdSphs(c, ■'^x) ~ 


(28) 


where (crn)**^ is the total annihilation cross section today as predicted by the model being 
tested with couplings c and DM mass m^. 


Monojets at LHC 

Here we take the most recent 95% CL LHC results (rescaled to the 68 % CL) from 
monojetanalyses [ 21 ], applied to the effective vector interaction operator qx^qxifiXy 
where a universal coupling to up and down-type quarks of both chiralities was assumed.® 

® In view of the possible issues regarding the validity of the EFT approach in collider searches [IS], we 
consider the limits from which respect the validity criteria of EFT. In the next LHC run, the issues 


14 





Notice that, if opposite couplings are assumed for the different chiralities, the same bound 
is recovered [SI]. Indeed, any interference effect between the different effective operators 
will be chirality suppressed. Thus, it is safe to assume that the total contribution can be 
obtained as the incoherent sum of the individual ones. Since a full collider simulation is 
beyond of the scope of this work, we recast the existing limits on the coefficient of 

the effective operator as a function of the DM mass. Since the constraint mainly stems 
from events with large Q^, for which the valence quarks dominate the parton distribution 
function of the proton, we assume that only the hrst generation quarks participate. Thus, 
the contribution to the function will be given by 

XlHc(C;''^x) = J'(^exp^)2^ ) (29) 

where = 3c^i + 2 c^r + 

Mono-photons at LEP 

Similarly to the monojet case, the 95% CL constraint (rescaled to the 68% CL) on 
the coefficient of the operator, by LEP mono-photon+ETgg searches is 

considered [8]. The contribution to the is thus: 

Xlep(c,’^x) ~ ’ (^^) 

where (C*'^p)2 = cL + c2^. 

IV. NUMERICAL DETAILS 

All the phenomenological models described in Sec. have been implemented in 
CalcHEP [52] using FeynRules [53]. The result is then fed into micrOmegas [35] for the 
computation of the DM observables. 

With the chi-square function implemented as described above, the parameter spaces of 
the different models need to be efficiently explored in order to derive constraints. While a 
simple grid scan may be viable in the case of the family-oriented model (with only three 

with the validity of the EFT will be even more relevant, and therefore an analysis in terms of Simplified 
Models [151150] might be preferable. 
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parameters), it quickly becomes prohibitively expensive for models with a larger number 
of free parameters. In particular, the general model with 15 free parameters would be 
unfeasible to scan in this fashion. Moreover, we expect the parameter space to be affected 
by several degeneracies. These further decrease the effectiveness of a grid scan and require 
the grid spacing to be very small in order to hnd the global minimum. For these reasons, 
we have opted to perform our scans using the nested sampling m Monte Carlo algorithm 
implemented in the MultiNest software [55]. This method is particularly designed for 
handling parameter space degeneracies as well as for preferentially scanning the regions of 
parameter space where the likelihood function C = exp(—x^/2) is large. Although the 
MultiNest software was initially designed for computing Bayesian evidence, it produces a 
sample of the parameter space with the corresponding likelihood values as a by-product. In 
this paper, we take a purely frequentist approach and only consider the values obtained 
from the likelihood. Thus, MultiNest is used only for its capability of efficiently sampling 
the regions with relatively low values. 

In our MultiNest simulations, the number of live points was set to 750, the tolerance 
to 0.2 and the efficiency to 0.9. Constant efficiency mode was not used. The range of 
parameters scanned was adapted to the values of the coefficients required to reproduce the 
correct relic abundance, for each of the dark matter masses considered. The scan was done 
in linear scale over both positive and negative values for the coefficients accompanying each 
of the operators. The number of distinct samples obtained ranges between 10® and 5 x 10®, 
depending on the particular model and DM mass considered. 


V. RESULTS 


This section describes the results obtained from a global £t using the models described 
in Sec. m All constraints considered in Sec. uni have been included in our simulations. 

Note that the dimensionful couplings Cj will depend on the effective theory mass scale A 
as A“^. In this section, we will use units of TeV“^ for the couplings Cj and the rough limit 
A^ < 1/cj may be used to assess the range of validity of the EFT assumption. 


Using the x^ function built as described in Sec. we have derived two main types of 
constraints on the parameters of the model: 


i) For a hxed DM mass m^, we scan the parameter space for all couplings in each 
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model. This results in a sample of points in the parameter space, which is generally 
concentrated in the regions corresponding to the lowest values of the function. We 
use these to perform parameter estimation and determine the allowed regions for the 
different couplings within a given model, for given values of m^. 

ii) The calculation described in (i) can be repeated for several values of the DM mass. 
By minimizing the global x^ over all couplings c in the model, the resulting Ay^, 

= min y2(c,mx) - miny2(c, m^), (31) 

C 

allows parameter estimation of the DM mass m^. By deriving conhdence intervals 
around the global minimum of a given model, this will indicate whether there are any 
values of the DM mass which are disfavoured by current data. Since the simulations 
are computationally rather expensive, they have been performed for a few values of 
the DM mass, namely = 10, 50, 100, 200, 500 and 1000 GeV. Our minimization 
in the variable is therefore performed in a discrete fashion, using only this set of 
values. However, from > 200 GeV on, we hnd essentially no change in the 
value. 


Our main results are summarized in Figs. and Figure shows the value of the 
Ay^ (see Eq. iHi) obtained for the different models under consideration, as indicated in 
the legend, see item (ii) above. The dashed horizontal line shows the value of the Ay^ 
corresponding to 3a limit for 1 d.o.f.. As can be seen from this hgure, DM masses around 
40 GeV and below are disfavored at 3cr with respect to higher masses for most models under 
consideration. The only exceptions are the general and leptophilic models, for which only 
masses below 20 GeV are disfavored. 

The fact that very light DM masses are generally disfavored can be understood from the 


complementarity between different data sets. As can be seen from Eq. (16), the annihilation 
cross section is proportional to m^, while the energy density scales as m^. Thus, models with 
very light DM masses will require large couplings to the SM in order to satisfy relic density 
constraints. However, since there is no positive signal from DM in collider, direct or indirect 
detection data, a signihcant tension between the different data sets occurs, which eventually 
increases the minimum value of the y^. The tension is stronger in models where DM either 
does not couple to leptons (leptophobic), or in models where the lepton couplings are related 
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(GeV) 


FIG. 1: Value of the Ax^ (see Eq. @) obtained for the different models under consideration (indicated 
in the legend) as a function of the DM mass in GeV. For reference, the value of the Ax^ corresponding to 
the 3cr limit (for f d.o.f.) is also shown by the dashed horizontal line. 


to others (i.e., flavour-blind or family-oriented). For the general and leptophilic models, the 
tension is relaxed since the bounds which mainly constrain the couplings to leptons (LEP, 
CMB and AMS data) are not as strong as those constraining the quark couplings (direct 
detection, LHC and FermiLAT bounds). Finally, one can observe a slight preference in the 
fit for DM masses around 100 GeV for all models except for the leptophobic and family- 
oriented. This can be understood from the mild preference of AMS data for a DM signal 
around this mass, since smaller masses are too strongly constrained, while heavier masses do 
not produce an appreciable signal. It should be stressed out, however, that this preference is 
far from being statistically significant. This feature is absent in the leptophobic and family- 
oriented models since the required fermion couplings for this signal to take place are either 
absent or constrained to be very small by other experimental bounds. 

The second main result of this paper is shown in Fig. where the allowed regions at 1, 2 
and 3a (for 1 d.o.f.) are shown for all couplings associated to the effective operators and for 
the five models under consideration, see Sec. For each coupling, our results are shown for 
three different values of the DM mass, 10 GeV (upper/blue bands), 100 GeV (middle/green 
bands) and 500 GeV (lower/red bands). Each panel also shows the constraint imposed by 
relic density on the weighted sum of the squares of all couplings, i.e.. Sc as defined in 


Eq. (16). 
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FIG. 2: Allowed absolute values for the couplings associated to the different operators present in each 
of the models considered in this work. Dark, medium and light shaded areas correspond to the allowed 
regions at 1, 2 and Scr, respectively, for 1 d.o.f.. For each coupling, the results shown by the upper/blue, 
middle/green and lower/red bands correspond to DM masses of 50, 100 and 500 GeV, respectively. The gray 
bands indicate that the coupling is unconstrained since the process XX ^ is not kinematically allowed for 
the considered value of the DM mass. 


The tension between different data sets can be nnderstood from the resnlts fonnd in Fig. 
As already explained, the constraint from relic density tends to bring the conplings to larger 
valnes as the DM mass is decreased. This is observed when comparing the upper/bine and 
lower/red bands, for all models under consideration, and in particular for the allowed values 
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of Sc. One can clearly see how the tension between different data sets can lead to having only 
one coupling as the major contributor to the relic density for a given model. For instance, 
for the general and leptophilic models with = 50 GeV, the relic density constraint is 
satished with a sizable coupling to the second family lepton doublet, and the preferred region 
for this coupling at la does not include zero. This can be understood as follows. First, the 
bounds on couplings to quarks of all generations are very strongly constrained by direct 
detection and Fermi-LAT experiments. Moreover, on the leptonic side, the bounds on RH 
leptons from both AMS (electrons) and Fermi-LAT (tans) are very strong. In addition, LH 
couplings would also imply DM annihilations through muon neutrinos, so that the indirect 
searches are relaxed while the relic abundance constraint is more easily fulhlled. As such, 
the weakest global bound appears for the second lepton generation, which can accommodate 
the required relic abundance without being in tension with AMS or Fermi-LAT. 

Nevertheless, a signihcant tension persists, as seen in Fig. This is particularly the case 
for light DM masses, see, e.g., the bands corresponding to 50 GeV in Fig. However, for 
100 GeV mass (and heavier) the situation is different. For such large values, the indirect 
detection constraints are already weaker than the relic abundance one. Thus, since a coupling 
to LH fermions implies two annihilation channels, both contributing to decrease the DM 
density, we see that the LH couplings are more constrained than the corresponding RH 
couplings. Similarly, due to color multiplicity factors, couplings to quarks are more strongly 
constrained than those to leptons (in particular the LH ones). 

Another interesting example is found in the panel for the flavour-blind model. Since 
the coupling to the muon is equal to that of the electron and tan (which are very strongly 
constrained) the relic abundance cannot be obtained out of annihilation to the second lepton 
doublet, as for the previous cases. For = 50 GeV, we see from the hgure that the DM 
would prefer to annihilate to up-type RH quarks instead (the Icr region for this coupling is 
around the value of 1). However, this results in a strong tension with the Fermi-LAT data, 
and this mass is therefore disfavored at ~ Sa with respect to the best £t at higher mass (see 
FigQ. 


Finally, as was already mentioned in Sec. |IH[ we have found two interesting degeneracies 
among different parameters in the models considered in this work. The hrst one is due to 
the constraints coming from direct detection experiments, which allows for a degeneracy 
between the couplings to quarks. This is shown in the left panel in Fig. where the allowed 
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FIG. 3: Left: Allowed regions in the parameter space corresponding to the couplings to the first family 
of quarks for the general model. The degeneracy between cqi,Cur and c^r is explicitly shown, see text 
for details. The weights of each coefficient are chosen so as to show explicitly the degeneracy condition in 


Eq. (18). Right: Allowed regions for the couplings to the second and third family in the family oriented 


scenario. The point C 2 = C 3 = 0 is clearly disfavored due to the efficient combination of relic abundance and 
direct detection constraints, see text for details. In both panels, the different colors correspond to different 
confidence levels, indicated in the legend, and the DM mass is = 100 GeV. 

confidence regions for the cqi, Cur and CdR conplings are depicted, showing the degeneracy 
explicitly. For = 100 GeV, the most stringent bonnd on the effective operators involving 
the hrst generation of quarks comes from direct detection experiments. Nevertheless, as can 
be seen from this figure, the allowed regions extend to rather large values of the couplings, 
as long as the relation 2.12cqi ~ —Cur — l.l2cdR is satisfied, given the values of Z and 


A of the Xenon material used by LUX (see Sec. Ill for details). While the inclusion of 
EDELWEISS-II, with a different target material, could in principle lift this degeneracy, the 
difference in the ratio of protons and neutrons is not large enough for this task. The hnal 
limits found in Fig. [^for the degeneracy line are rather stemming from the relic abundance 
constraint. 

The second degeneracy stems from the strong constraints from relic abundance on Sc, 
which generally imply that all couplings must lie on the surface of a hyperellipsoid. In 
particular, for the family oriented scenario, this reduces to the ellipse displayed in the right 
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panel in Fig. where the allowed regions at 1, 2 and 3a are shown in the C 2 — C 3 plane for 
2 d.o.f.. Given that the strongest constraints from colliders and direct detection experiments 
apply to the hrst generation, and in this model all particles in the hrst generation have the 
same coupling, the constraints on ci are very strong since the direct detection degeneracy 
relation cannot be fulhlled. This is shown in Fig. where it can be seen that the coupling 
to the hrst family is restricted to ci < 0.001/TeV^ at 2 ct for this model. This model only has 
two additional couplings, C 2 and C 3 . Therefore, in order to satisfy relic density constraints, 
either the coupling to the second or the third family (or both) have to be different from 
zero. As a consequence, the allowed region is shaped as an ellipse in the C 2 — C 3 plane. 

A hnal remark regarding the validity of the EFT is in order. As can be seen from 
eq. (16) and from Fig. the weighted sum of the coefficients required to obtain the correct 
thermal relic abundance goes like Sc oc l/rriy.. On the other hand, for DM annihilation, 
the EFT validity condition reads A > 2m^. Since Sc is a combination of couplings, we can 
generically write it as Sc ~ 1/A^. Therefore, by combining the validity condition with the 
relic abundance requirement, a maximum DM mass is obtained for which the EFT ceases 
to be valid about 0{2 — 5 TeV) . Notice however that we do not expect the data to strongly 
constrain thermal DM in this regime. 


VI. CONCLUSIONS 

In this work we have explored how strongly the combination of present data from very 
different experimental probes can constrain the different couplings between Dark Matter 
(DM) and the Standard Model (SM) particle content. The focus of the work was a bottom- 
up approach to understand what data itself is able to tell us regarding DM interactions, 
minimizing when possible any theoretical input or bias. To this end, we make use of an 
Effective Field Theory (EFT) approach to attain the desired model-independence at the 
expense of assuming that DM-SM interactions are mediated by heavy particles that can be 
integrated out of the theory. Furthermore, we have allowed independent couplings of DM to 
all SM helds so as to minimize any theoretical bias, but we have limited these interactions to 
flavour-conserving ones in the SM sector, since generally stronger constraints apply to them. 
Finally, if all possible Lorentz structures were simultaneously allowed, the parameter space 
would become too large to explore and constrain with present data. Thus, we restrict our 


22 



analysis to Dirac DM fermions which interact with the SM constituents through independent 
operators of the form 

Ci,p ixYx) , (32) 

with independent coefficients Cj^p for all SM fermions /j and chiralities P = Pp, Pp. 

This working model, which is not intended to reproduce any particular ultraviolet comple¬ 
tion, provides the necessary parametrization to assess through a joint analysis the constraints 
on the individual interaction strengths from present data. We constrain 15 independent co¬ 
efficients through a combination of 8 experimental probes comprising relic abundance, direct 
and indirect DM searches as well as collider limits. In addition to this general setup, we 
also explore how these constraints are affected when correlations between the different coef¬ 
ficients are allowed, or when some of the operators are forbidden. In particular, we studied 
the completely leptophilic, leptophobic and flavour-blind cases, as well a family-oriented 
scenario in which all members of the same generation share the same coupling to the DM 
particle. 

From our global fits we find that for DM masses < 20 GeV the tension with respect to 
the best fit between the couplings necessary to reproduce the observed DM relic abundance 
and the upper bounds from null results exceeds the 3a level, for the general and leptophilic 
scenarios considered, while the same happens for < 40 GeV for the leptophobic, flavour 
blind and family oriented cases. For these three cases, the rises very steeply as the DM 
mass decreases such that for < (9(10 — 20) GeV this tension reaches the 5a level, within 
our EFT framework. 

Furthermore we find that, due to the slightly weaker present constraints, couplings to 
the second-generation lepton doublet are preferred, whenever they are available within the 
model under study. If this coupling is not available or is instead further constrained by an 
assumed correlation, as is the case of the flavour-blind scenario, a more sizable coupling 
to the right-handed, up-type quark singlet is instead preferred for similar reasons. Thus, 
it would be interesting to improve our present constraints on these couplings. Finally, we 
also find that, since all couplings are assumed to be independent and free to vary in the fit, 
the very stringent constraints stemming from direct search experiments such as LUX imply 
that either the first generation quark couplings to DM are extremely small or that they are 
related to each other such that the different contributions to these processes cancel against 
each other, i.e., 2 cqi -f- Cur + CdR ~ 0. 
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DM searches worldwide are now probing and constraining essentially all possible interac¬ 
tion channels between DM and the known matter constituents through extremely different 
and complementary and search techniques. In this context, and given our lack of a unique 
theoretical DM paradigm, it is important to test different DM models against present data 
with bottom-up approaches where theoretical biases are not imposed. Unfortunately, true 
and complete model independence cannot be achieved. Indeed, while EFT offers an ideal 
frame for this sort of studies, its adoption already enforces some assumptions, such as the de¬ 
coupling nature of the mediating particle or the uniqueness of the DM candidate. Moreover, 
true model independence would imply the inclusion of hundreds of independent operators 
with different flavour and chiral structures, rendering the analysis too general to actually 
provide any useful information. In this work we have reduced the number of theoretical as¬ 
sumptions taking a hrst step which implies some unavoidable restrictions to the number and 
types of effective operators included. It would be very interesting to supplement our results 
with additional complementary analyses, including different Lorentz structures, different 
DM helds (Majorana fermions, or scalars), or with extensions beyond the EFT approach by 
allowing one (or several) generic light mediator(s). 
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Appendix A: Fitting procedure for the AMS02 positron flux data 


As already explained in Sec. the measurements from AMS on the positron fraction 
are also used to derive upper bounds on DM annihilation cross sections. The main issue that 
has to be dealt with when doing so, however, are the large uncertainties on the positron and 
electron fluxes from astrophysical sources. It is common to use a linear combination of two 
power laws to parametrize them (see, e.g.. Refs |12l ESI EZ]), where typically the positron 
flux includes an exponential cut-off at high energies. Our parametrization for the electron 


and positron fluxes is shown in Fq. (23). For the propagation of charged cosmic rays in the 


astrophysical media, we use the MED parameters dehned in Ref. 

In principle, both DM and the astrophysical contribution to the electron and positron 
fluxes are affected by the solar modulation. This can be explicitly taken into account by 
computing the flux at the top of the atmosphere (©) under the force field approximation. 


see Eq. (24). Therefore, under these assumptions our background already depends on 9 


(Al) 


parameters 

{Ce,7e,C^,7s,F^,(/>+} for positrons, and 

{Ci, 7 i,C 2 , 72 , 0 “} for electrons. 

To get the total positron fraction, one would add to the background the additional contribu¬ 
tion from DM annihilation, which in principle depends on the solar modulation parameter 


as well, see Eq. (25) 


In principle, the best thing would be to perform a fit to the AMS02 data where, for 
each value of the DM mass, the minimum of the is searched for after marginalizing 
over the 9 parameters listed above. However, this is computationally rather expensive. 
Furthermore, we found that severe degeneracies take place among the different parameters, 
which complicates the problem even further. 

However, the problem can be considerably simplified by considering the following. In 
principle, the positron flux will be the most sensitive to the DM annihilation signal, while 
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the electron flux will be mainly dominated by astrophysical backgrounds instead. This allows 
to reduce the number of parameters in the fit significantly: since the electron flux will be 
independent from the signal, it can be fitted independently and leave the parameters fixed 


during the fit. This is done using the parametrization in Eq. (23) for the electron flux, and 


the publicly available data® on the electron flux from the Fermi LAT collaboration 

When fitting the Fermi LAT data, a fit is performed considering both the low-energy 
and high-energy data sets, between 7 GeV and 1 TeV. The resulting curve is shown in the left 
panel in Fig. together with the data points as extracted from Ref. [ 13 ]. The uncertainties 
in each bin are computed by adding in quadrature their statistical and systematic errors.^® 
We And that the parameters which give a best-fit to the Fermi LAT electron flux data are 


^1(2) = 213.7(140.1) s-Wr-i m-2 GeV^ 

71(2) = 3.7 (3.0); r =0.0. 

For these values of the parameters, we And a minimum y^/d.o.f. = 4.7/38. 


(A2) 



FIG. 4: (Left) Our best-fit for the electron flux compared to Fermi LAT data. (Right) Our best-fit to the 
positron fraction in the absence of any contribution from DM annihilation, compared to the AMS data. See 
text for details. 

Finally, when deriving the constraint from AMS data, the total positron flux is obtained 
as the addition of the astrophysical background (which depends on Ge, Gs, 7 e, 7 s, and 

® The AMS02 data in Ref. [TS] contains only the positron fraction and fluxes, but not the electron data. 
Therefore, we use the Fermi LAT data in this case, which is publicly available |43) . 

In the case of asymmetric systematic errors, we (conservatively) take the largest value. 
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0+) and the contribution from DM annihilation. In the absence of an extra contribution 
from DM annihilation, we hnd that the following parameters give a best-£t to the positron 
fraction data from AMS 


Ce{s) = 30.4(2.0) s ^sr ^m ^ GeV 

7e(s) = 3.9 (2.5); = 1086.8GeV; 0+ = 0.0, 


(A3) 


with x^/d.o.f. = 26.2/53. The positron fraction obtained with these parameters can be seen 
in the right panel in Fig. [^together with the AMS data. Again in this case, the errors are 
taken as the sum in quadrature of the statistical and systematic errors in each bin. 

In our simulations, however, we include the DM annihilation to the positron flux and 


we let the parameters in Eq. (Al) vary during the £t. When doing so, we hnd that the 
solar modulation parameter does not have a major impact in the fit while it generates some 
numerical degeneracies with other parameters, which are difficult to deal with. Therefore, 
we also £x this parameter to the value which gives a best-£t to the AMS02 data using the 


background contribution alone. The rest of the parameters in Eq. (Al) are left free during 
the £t, and will be htted for each value of the DM mass and the couplings in an independent 
way. 



FIG. 5: Limits to the DM interaction cross section as a function of the DM mass for several primary 
annihilation channels to leptons, as indicated in the legend, when only one coupling is allowed at a time. 
The regions above each line are excluded at 90% CL (1 dof). The horizontal line indicates the cross section 
needed to satisfy relic abundance constraints. 

Finally, even though in our simulations we consider several couplings between the DM 
and SM particles at once, it is useful to look at the limits obtained when only one coupling 
to a SM fermion is allowed at a time. This allows to illustrate the interplay between different 
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FIG. 6; Same as Fig. [^but for primary DM annihilation into quarks, as indicated in the legend. 

data sets and where the tension in the £t for low values of the DM mass may come from. 
The limits on the DM-SM interaction cross section are shown in Fig. for leptons and in 
Fig. for quarks, as a function of the DM mass. Since in our £t we are combining AMS 
with relic abundance constraints, a signihcant tension will only take place when the limit on 
the interaction cross section gets below the value needed for relic abundance. Therefore, one 
can already see from this figure that the AMS data will be most effective in constraining the 
couplings of the DM to electrons (for < 10 GeV) and muons (for < 60 GeV), but 
will be less efficient for other DM fermions (for instance, for tans it will only significantly 
affect the £t for < 20 GeV). 
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